# ================================================================================== #
# Complaints About Police Misconduct Have Adverse Effects for Black Civilians (PSRM)
# - Script: Main analyses reported in manuscript
# - Author: Patrick Kraft (patrickwilli.kraft@uc3m.es)
# ================================================================================== #


# Load raw data, packages, and custom functions ---------------------------

source(here::here("00-func.R"))
load(here("sqf_ccrb.Rdata"))


# Figure 1: Overview of time series ---------------------------------------

p1a <- ggplot(total, aes(x=date, y=sqf)) + plot_default + 
  geom_line() + ylab("Number of Incidents (x 1000)") + xlab("Time") +
  ggtitle("(A) Monthly Stop, Question, Frisk Incidents")

p1b <- ggplot(total, aes(x=date, y=ccrb)) + plot_default +
  geom_line() + ylab("Number of Complaints (x 100)") + xlab("Time") +
  ggtitle("(B) Monthly CCRB Complaints")

ggsave("out/fig1_overview.png", 
       grid.arrange(p1a,p1b,ncol=1), 
       height=4.5, width=6.5)


# Figure 2: VAR effects of CCRB complaints on subsequent SQF --------------

m2 <- NULL

## variable selection for VAR
endo <- c("sqf_d","ccrb_d","arrests_d")
exo <- c("media_d", "unemp_d","visit_overseas_d",
         "temp_mean","precip_total", "jan")

## all incidents
tmp <-  na.omit(left_join(total, controls)[, c(endo, exo)])
m2$total <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")

## white complainants / suspects
tmp <-  na.omit(left_join(white, controls)[, c(endo, exo)])
m2$white <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")

## black complainants / suspects
tmp <-  na.omit(left_join(black, controls)[, c(endo, exo)])
m2$black <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")

## hispanic complainants / suspects
tmp <-  na.omit(left_join(hispanic, controls)[, c(endo, exo)])
m2$hispanic <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")

## prepare estimates for plot
m2_coef <- map_dfr(m2, extractCoefs, .id = "id") %>%
  mutate(p = round(2* pt(abs(Estimate/SE), df=m2[[1]]$obs, lower.tail = F), 3),
         CCRB = factor(gsub(".*\\.l","Lag ", Variables), levels=paste0("Lag ",4:1)),
         Group = factor(id, levels = c("total", "white", "black", "hispanic"), 
                        labels = c("All CCRB Complaints / SQF Incidents","Whites (in CCRB & SQF)",
                                   "Blacks (in CCRB & SQF)","Latinos (in CCRB & SQF)")))

## coefficient plot
p2a <- ggplot(m2_coef, aes(y=Estimate, ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE, x=CCRB)) + 
  plot_default + geom_linerange() + geom_point(size=1) + coord_flip() + 
  facet_wrap(~Group, ncol=1, scales="free_y") +
  geom_text(label = paste0("p = ", m2_coef$p), nudge_x = .4, size = 2) +
  geom_hline(yintercept=0, linetype="dotted") + xlab("CCRB Complaints") +
  ggtitle("(A) Estimated Effect on SQF Incidents") + ylab("Coefficient")

## simulate series
m2_series <- rbind(
  varSim(m2[[1]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)), 
         shock = list(ccrb_d = max(total[,"ccrb_d"],na.rm = T)), 
         level = apply(total[,c("sqf","ccrb","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "All CCRB Complaints / SQF Incidents", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group), 
  varSim(m2[[2]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_d = max(white[,"ccrb_d"],na.rm = T)),
         level = apply(white[,c("sqf","ccrb","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Whites (in CCRB & SQF)", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group),
  varSim(m2[[3]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_d = max(black[,"ccrb_d"],na.rm = T)),
         level = apply(black[,c("sqf","ccrb","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Blacks (in CCRB & SQF)", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group),
  varSim(m2[[4]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_d = max(hispanic[,"ccrb_d"],na.rm = T)),
         level = apply(hispanic[,c("sqf","ccrb","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Latinos (in CCRB & SQF)", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group)) %>%
  mutate(Group = factor(Group, levels = c("All CCRB Complaints / SQF Incidents",
                                          "Whites (in CCRB & SQF)",
                                          "Blacks (in CCRB & SQF)",
                                          "Latinos (in CCRB & SQF)")))

## compute confidence intervals for simulations
m2_ci <- m2_series %>% group_by(time, Group) %>%
  summarize(avg=mean(SQF), cilo=quantile(SQF,.025), cihi=quantile(SQF,.975), SQF=min(SQF))
m2_text <- m2_ci %>% 
  group_by(Group) %>% summarize(SQF=min(SQF)) %>%
  mutate(shock = c(max(total[,"ccrb_d"],na.rm = T), max(white[,"ccrb_d"],na.rm = T),
                   max(black[,"ccrb_d"],na.rm = T), max(hispanic[,"ccrb_d"],na.rm = T)),
         shock = paste0(" + ",round(shock, 1) * 100," complaints"),
         time = 0)

## plot simulations
p2b <- ggplot(m2_series, aes(x=as.factor(time), y=SQF)) + plot_default +
  geom_line(aes(group=iter), alpha=.01) + ylab("Number of Incidents (x 1000)") + 
  geom_vline(xintercept = 2) + facet_wrap(~Group,ncol=1,scales="free_y") + xlab("Month") + 
  ggtitle("(B) Predicted Monthly SQF Incidents") +
  geom_linerange(aes(ymin=cilo, ymax=cihi, x=as.factor(time)), data=m2_ci[m2_ci$time>0,]) +
  geom_point(aes(y=avg, x=as.factor(time)), size=1, data=m2_ci[m2_ci$time>0,]) +
  geom_segment(aes(x = 2, y = avg, xend = 6, yend = avg), lty="dotted", data=m2_ci[m2_ci$time==0,]) +
  geom_text(aes(label=shock), vjust=0, hjust=0, size=2, data=m2_text)

## save combined plot
ggsave("out/fig2_main.png", 
       grid.arrange(p2a, p2b, ncol=2), 
       height=5.5, width=6.5)


# Figure 3: VAR effects by officer race (black suspects/complainants) --------

m3 <- NULL

## variable selection for VAR
exo <- c("media_d", "unemp_d","visit_overseas_d",
         "temp_mean","precip_total", "jan")

## by white/other police officer
series <- c("po","pw")
for(i in series){
  endo <- c("sqf_d",paste0("ccrb_",i,"_d"),"arrests_d")
  tmp <- na.omit(left_join(black, controls)[, c(endo, exo)])
  m3[[i]] <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")
}

## prepare estimates for plot
m3_coef <- m3 %>% 
  map_dfr(~filter(prePlot(.$varresult, NeweyWest = F, report = "ccrb"), Model=="sqf_d")) %>%
  mutate(p = round(2* pt(abs(Estimate/SE), df=m3[[1]]$obs, lower.tail = F), 3),
         CCRB = factor(gsub(".*\\.l","Lag ", Variables), levels=paste0("Lag ",4:1)),
         Group = factor(gsub("_d\\..*","", Variables), 
                        levels = c("ccrb_po","ccrb_pw"), 
                        labels = c("Complaints against Non-White Officers",
                                   "Complaints against White Officers")))

## coefficient plot
p3a <- ggplot(m3_coef, aes(y=Estimate, ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE, x=CCRB)) + 
  plot_default + geom_linerange() + geom_point(size=1) + coord_flip() + 
  facet_wrap(~Group, ncol=1, scales="free_y") +
  geom_text(label = paste0("p = ", m3_coef$p), nudge_x = .4, size = 2) +
  geom_hline(yintercept=0, linetype="dotted") + xlab("CCRB Complaints by Blacks") +
  ggtitle("(A) Estimated Effect on SQF Incidents") + ylab("Coefficient")

## simulate series
m3_series <- rbind(
  varSim(m3[[1]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_po_d = max(black[,"ccrb_po_d"],na.rm = T)),
         level = apply(black[,c("sqf","ccrb_po","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Complaints against Non-White Officers", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group),
  varSim(m3[[2]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_pw_d = max(black[,"ccrb_pw_d"],na.rm = T)),
         level = apply(black[,c("sqf","ccrb_pw","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Complaints against White Officers", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group)) %>%
  mutate(Group = factor(Group, levels = c("Complaints against Non-White Officers", 
                                          "Complaints against White Officers")))

## compute confidence intervals for simulations
m3_ci <- m3_series %>% group_by(time, Group) %>%
  summarize(avg=mean(SQF), cilo=quantile(SQF,.025), cihi=quantile(SQF,.975), SQF=min(SQF))
m3_text <- m3_ci %>% 
  group_by(Group) %>% summarize(SQF=min(SQF)) %>%
  mutate(shock = c(max(black[,"ccrb_po_d"],na.rm = T), max(black[,"ccrb_pw_d"],na.rm = T)),
         shock = paste0(" + ",round(shock, 1) * 100," complaints"),
         time = 0)

## plot simulations
p3b <- ggplot(m3_series, aes(x=as.factor(time), y=SQF)) + plot_default +
  geom_line(aes(group=iter), alpha=.01) + ylab("Number of Incidents (x 1000)") + 
  geom_vline(xintercept = 2) + facet_wrap(~Group,ncol=1) + xlab("Month") + 
  ggtitle("(B) Predicted Monthly SQF Incidents") + 
  geom_linerange(aes(ymin=cilo, ymax=cihi, x=as.factor(time)), data=m3_ci[m3_ci$time>0,]) +
  geom_point(aes(y=avg, x=as.factor(time)), size=1, data=m3_ci[m3_ci$time>0,]) +
  geom_segment(aes(x = 2, y = avg, xend = 6, yend = avg), lty="dotted", data=m3_ci[m3_ci$time==0,]) +
  geom_text(aes(label=shock), vjust=0, hjust=0, size=2, data=m3_text)

## save combined plot
ggsave("out/fig3_orace.png", 
       grid.arrange(p3a, p3b, ncol=2), 
       height=3.5, width=6.5)


# Figure 4: Officer knowledge & case disposition (black suspects/complainants) --------

m4 <- NULL

## by substantiated / not substantiated
series <- c("subst1","subst2","subst3")
for(i in series){
  endo <- c("sqf_d",paste0("ccrb_",i,"_d"),"arrests_d")
  tmp <- na.omit(left_join(black, controls)[, c(endo, exo)])
  m4[[i]] <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")
}

## prepare estimates for plot
m4_coef <- m4 %>% 
  map_dfr(~filter(prePlot(.$varresult, NeweyWest = F, report = "ccrb"), Model=="sqf_d")) %>%
  mutate(p = round(2* pt(abs(Estimate/SE), df=m4[[1]]$obs, lower.tail = F), 3),
         CCRB = factor(gsub(".*\\.l","Lag ", Variables), levels=paste0("Lag ",4:1)),
         Group = factor(gsub("_d\\..*","", Variables), 
                        levels = c("ccrb_subst1","ccrb_subst2","ccrb_subst3"), 
                        labels = c("Officer not Contacted by CCRB",
                                   "Officer Contacted & Unsubstantiated by CCRB",
                                   "Officer Contacted & Substantiated by CCRB")))

## coefficient plot
p4a <- ggplot(m4_coef, aes(y=Estimate, ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE, x=CCRB)) + 
  plot_default + geom_linerange() + geom_point(size=1) + coord_flip() + facet_wrap(~Group, ncol=1) +
  geom_text(label = paste0("p = ", m4_coef$p), nudge_x = .4, size = 2) +
  geom_hline(yintercept=0, linetype="dotted") + xlab("CCRB Complaints by Blacks") +
  ggtitle("(A) Estimated Effect on SQF Incidents") + ylab("Coefficient")

## simulate series
m4_series <- rbind(
  varSim(m4[[1]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_subst1_d = max(black[,"ccrb_subst1_d"],na.rm = T)),
         level = apply(black[,c("sqf","ccrb_subst1","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Officer not Contacted by CCRB", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group), 
  varSim(m4[[2]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_subst2_d = max(black[,"ccrb_subst2_d"],na.rm = T)), 
         level = apply(black[,c("sqf","ccrb_subst2","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Officer Contacted & Unsubstantiated by CCRB", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group), 
  varSim(m4[[3]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_subst3_d = max(black[,"ccrb_subst3_d"],na.rm = T)),
         level = apply(black[,c("sqf","ccrb_subst3","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Officer Contacted & Substantiated by CCRB", SQF = sqf_d.1) %>%
    dplyr::select(SQF,time,iter,Group)) %>%
  mutate(Group = factor(Group, levels = c("Officer not Contacted by CCRB",
                                          "Officer Contacted & Unsubstantiated by CCRB",
                                          "Officer Contacted & Substantiated by CCRB")))

## compute confidence intervals for simulations
m4_ci <- m4_series %>% group_by(time, Group) %>%
  summarize(avg=mean(SQF), cilo=quantile(SQF,.025), cihi=quantile(SQF,.975), SQF=min(SQF))
m4_text <- m4_ci %>% 
  group_by(Group) %>% summarize(SQF=min(SQF)) %>%
  mutate(shock = c(max(black[,"ccrb_subst1_d"],na.rm = T), 
                   max(black[,"ccrb_subst2_d"],na.rm = T),
                   max(black[,"ccrb_subst3_d"],na.rm = T)),
         shock = paste0(" + ",round(shock, 1) * 100," complaints"),
         time = 0)

## plot simulations
p4b <- ggplot(m4_series, aes(x=as.factor(time), y=SQF)) + plot_default +
  geom_line(aes(group=iter), alpha=.01) + ylab("Number of Incidents (x 1000)") +
  geom_vline(xintercept = 2) + facet_wrap(~Group,ncol=1) + xlab("Month") +
  ggtitle("(B) Predicted Monthly SQF Incidents") +
  geom_linerange(aes(ymin=cilo, ymax=cihi, x=as.factor(time)), data=m4_ci[m4_ci$time>0,]) +
  geom_point(aes(y=avg, x=as.factor(time)), size=1, data=m4_ci[m4_ci$time>0,]) +
  geom_segment(aes(x = 2, y = avg, xend = 6, yend = avg), lty="dotted", data=m4_ci[m4_ci$time==0,]) +
  geom_text(aes(label=shock), vjust=0, hjust=0, size=2, data=m4_text)

## save combined plot
ggsave("out/fig4_subst.png", 
       grid.arrange(p4a, p4b, ncol=2), 
       height=4.5, width=6.5)


# Figure 5: Officer discretion for stop (black suspects/complainants) --------

m5 <- NULL

## by officer discretion
series <- c("invest0","invest1",
            "report0","report1",
            "proxim0","proxim1",
            "vcrime_cs0","vcrime_cs1",
            "sights0","sights1")
for(i in series){
  endo <- c(paste0("sqf_",i,"_d"),"ccrb_d","arrests_d")
  tmp <- na.omit(left_join(black, controls)[,c(endo, exo)])
  m5[[i]] <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")
}

## prepare estimates for plot
m5_coef <- map2_dfr(.x = m5, .y = paste0("sqf_",series,"_d"), 
                    ~filter(prePlot(.x$varresult, NeweyWest = F, report = "ccrb"), Model==.y)) %>%
  mutate(p = round(2* pt(abs(Estimate/SE), df=m5[[1]]$obs, lower.tail = F), 3),
         CCRB = factor(gsub(".*\\.l","Lag ", Variables), levels=paste0("Lag ",5:1)),
         Group = factor(Model, levels = paste0("sqf_",series,"_d"), 
                        labels = c("SQF due to Ongoing Investigation: No", 
                                   "SQF due to Ongoing Investigation: Yes",
                                   "SQF due to Report by Victim/Witness: No",
                                   "SQF due to Report by Victim/Witness: Yes",
                                   "SQF due to Proximity to Scene of Crime: No",
                                   "SQF due to Proximity to Scene of Crime: Yes",
                                   "SQF due to Observed Violent Crime: No",
                                   "SQF due to Observed Violent Crime: Yes",
                                   "SQF due to Sights/Sounds of Criminal Activity: No",
                                   "SQF due to Sights/Sounds of Criminal Activity: Yes")))

## coefficient plot
p5 <- ggplot(m5_coef, aes(y=Estimate, ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE, x=CCRB)) + 
  plot_default + geom_linerange() + geom_point(size=1) + coord_flip() + facet_wrap(~Group, ncol=2) +
  geom_text(label = paste0("p = ", m5_coef$p), nudge_x = .4, size = 2) +
  geom_hline(yintercept=0, linetype="dotted") + xlab("CCRB Complaints by Blacks") +
  ggtitle("Estimated Effect on SQF Incidents (Black Suspects)") + ylab("Coefficient")

## save plot
ggsave("out/fig5_discr.png", 
       p5, 
       height=6.5, width=6.5)


# Figure 6: Results by Borough (Black suspects only) ----------------------

m6 <- NULL

## By borough
boro <- c("Bronx","Brooklyn","Queens","Manhattan","Staten_Island")
boroLab <- c("Bronx","Brooklyn","Queens","Manhattan","Staten Island")
for(i in boro){
  endo <- c(paste0(c("sqf_","ccrb_"),i,"_d"),"arrests_d")
  tmp <- na.omit(left_join(black, controls)[,c(endo, exo)])
  m6[[i]] <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")
}

## prepare estimates for plot
m6_coef <- map2_dfr(.x = m6, .y = paste0("sqf_", boro, "_d"), 
                    ~filter(prePlot(.x$varresult, NeweyWest = F, report = "ccrb"), Model==.y)) %>%
  mutate(p = round(2* pt(abs(Estimate/SE), df=m6[[1]]$obs, lower.tail = F), 3),
         CCRB = factor(gsub(".*\\.l","Lag ", Variables), levels=paste0("Lag ",5:1)),
         Group = factor(Model, levels = paste0("sqf_",boro,"_d"), 
                        labels = paste0(boroLab," (CCRB & SQF)")))

## coefficient plot
p6a <- ggplot(m6_coef, aes(y=Estimate, ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE, x=CCRB)) + 
  plot_default + geom_linerange() + geom_point(size=1) + coord_flip() + 
  facet_wrap(~Group, ncol=1, scales="free_y") +
  geom_text(label = paste0("p = ", m6_coef$p), nudge_x = .4, size = 2) +
  geom_hline(yintercept=0, linetype="dotted") + xlab("CCRB Complaints") +
  ggtitle("(A) Estimated Effect on SQF Incidents") + ylab("Coefficient")

## simulate series
m6_series <- rbind(
  varSim(m6[[1]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_Bronx_d = max(black[,"ccrb_Bronx_d"],na.rm = T)), 
         level = apply(black[,c("sqf_Bronx","ccrb_Bronx","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Bronx (CCRB & SQF)", SQF = sqf_Bronx_d.1) %>%
    dplyr::select(SQF,time,iter,Group), 
  varSim(m6[[2]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_Brooklyn_d = max(black[,"ccrb_Brooklyn_d"],na.rm = T)), 
         level = apply(black[,c("sqf_Brooklyn","ccrb_Brooklyn","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Brooklyn (CCRB & SQF)", SQF = sqf_Brooklyn_d.1) %>%
    dplyr::select(SQF,time,iter,Group),
  varSim(m6[[3]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_Queens_d = max(black[,"ccrb_Queens_d"],na.rm = T)),
         level = apply(black[,c("sqf_Queens","ccrb_Queens","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Queens (CCRB & SQF)", SQF = sqf_Queens_d.1) %>%
    dplyr::select(SQF,time,iter,Group),
  varSim(m6[[4]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_Manhattan_d = max(black[,"ccrb_Manhattan_d"],na.rm = T)),
         level = apply(black[,c("sqf_Manhattan","ccrb_Manhattan","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Manhattan (CCRB & SQF)", SQF = sqf_Manhattan_d.1) %>%
    dplyr::select(SQF,time,iter,Group),
  varSim(m6[[5]], dumvar = apply(controls[,exo],2,function(x) rep(mean(x, na.rm=T), 6)),
         shock = list(ccrb_Staten_Island_d = max(black[,"ccrb_Staten_Island_d"],na.rm = T)),
         level = apply(black[,c("sqf_Staten_Island","ccrb_Staten_Island","arrests")],2,mean,na.rm=T)) %>%
    mutate(Group = "Staten Island (CCRB & SQF)", SQF = sqf_Staten_Island_d.1) %>%
    dplyr::select(SQF,time,iter,Group)) %>%
  mutate(Group = factor(Group, levels = paste0(boroLab," (CCRB & SQF)")))

## compute confidence intervals for simulations
m6_ci <- m6_series %>% group_by(time, Group) %>%
  summarize(avg=mean(SQF), cilo=quantile(SQF,.025), cihi=quantile(SQF,.975), SQF=min(SQF))
m6_text <- m6_ci %>% 
  group_by(Group) %>% summarize(SQF=min(SQF)) %>%
  mutate(shock = c(max(black[,"ccrb_Bronx_d"],na.rm = T), 
                   max(black[,"ccrb_Brooklyn_d"],na.rm = T), 
                   max(black[,"ccrb_Queens_d"],na.rm = T), 
                   max(black[,"ccrb_Manhattan_d"],na.rm = T), 
                   max(black[,"ccrb_Staten_Island_d"],na.rm = T)), 
         time = 0) %>% mutate(shock = paste0(" + ",round(shock, 1) * 100," complaints"))

p6b <- ggplot(m6_series, aes(x=as.factor(time), y=SQF)) + plot_default +
  geom_line(aes(group=iter), alpha=.01) + ylab("Number of Incidents (x 1000)") + geom_vline(xintercept = 2) +
  facet_wrap(~Group,ncol=1,scales="free_y") + xlab("Month") +
  ggtitle("(B) Predicted Monthly SQF Incidents") +
  geom_linerange(aes(ymin=cilo, ymax=cihi, x=as.factor(time)), data=m6_ci[m6_ci$time>0,]) +
  geom_point(aes(y=avg, x=as.factor(time)), size=1, data=m6_ci[m6_ci$time>0,]) +
  geom_segment(aes(x = 2, y = avg, xend = 6, yend = avg), lty="dotted", data=m6_ci[m6_ci$time==0,]) +
  geom_text(aes(label=shock), vjust=0, hjust=0, size=2, data=m6_text)

## save combined plot
ggsave("out/fig6_boro.png", 
       grid.arrange(p6a,p6b,ncol=2), 
       height=6.5, width=6.5)


# Figure 7: Results by precinct Black population size ---------------------

m7 <- NULL

## by officer discretion
series <- c("pctblack_Bronx", "pctblack_Brooklyn", "pctblack",
            "pctother_Bronx", "pctother_Brooklyn", "pctother")
for(i in series){
  endo <- c(paste0("sqf_",i,"_d"),paste0("ccrb_",i,"_d"),"arrests_d")
  tmp <- na.omit(left_join(black, controls)[,c(endo, exo)])
  m7[[i]] <- VAR(tmp[,endo], exogen = tmp[,exo], lag.max=5, ic="AIC")
}

## prepare estimates for plot
m7_coef <- map2_dfr(.x = m7, .y = paste0("sqf_",series,"_d"), 
                    ~filter(prePlot(.x$varresult, NeweyWest = F, report = "ccrb"), Model==.y)) %>%
  mutate(p = round(2* pt(abs(Estimate/SE), df=m7[[1]]$obs, lower.tail = F), 2),
         CCRB = factor(gsub(".*\\.l","Lag ", Variables), levels=paste0("Lag ",5:1)),
         Group = factor(Model, levels = paste0("sqf_",series,"_d"), 
                        labels = c("Percent Black Civilians\nAbove Median<>Bronx",
                                   "Percent Black Civilians\nAbove Median<>Brooklyn",
                                   "Percent Black Civilians\nAbove Median<>New York City",
                                   "Percent Black Civilians\nBelow Median<>Bronx",
                                   "Percent Black Civilians\nBelow Median<>Brooklyn",
                                   "Percent Black Civilians\nBelow Median<>New York City"))) %>% 
  separate(Group, into = c("percent","boro"), sep = "<>")

## coefficient plot
p7 <- ggplot(m7_coef, aes(y=Estimate, ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE, x=CCRB)) + 
  plot_default + geom_linerange() + geom_point(size=1) + coord_flip() + facet_grid(percent~boro, scales = "free_y") +
  geom_text(label = paste0("p = ", m7_coef$p), nudge_x = .4, size = 2) +
  geom_hline(yintercept=0, linetype="dotted") + 
  labs(y = "Lag", x = "Coefficient")

## save plot
ggsave("out/fig7_pctrace.png", 
       p7, 
       height=3, width=7)


# Save model objects for appendix ------------------------------------------------------

save(m2, m3, m4, m5, m6, m7, file = "out/results.Rdata")
